Packages
library(tidyverse) #Data manipulation and formatting
library(knitr) #Report tools and formatting
library(corrplot) #Pretty correlation plots
library(caret) #Multiuse package containing many models
library(sem) #Structural equation models for CFA
library(psych) #Exploratory factor models
library(randomForest) #Classification trees also in caret
library(gbm) #Logistic classifies for caret and other packages
library(DiagrammeR)
library(rattle)
library(Hmisc)
library(lavaan)
library(MVN)
library(GPArotation)
Data Import
BehDat.cntrl <- read.csv("Vole_Beh_R_Supp.csv") %>%
filter(TRT != "5MT") %>%
mutate(SP_INDEX_inv = SP_INDEX/-1) %>%
mutate(NO_TOT_DIST = NO_TRIAL2_TOT_DIST + NO_TRIAL3_TOT_DIST) %>%
mutate(bidose = ifelse(TRT=="CTRL", "CTRL","TRT")) %>%
select(
VOLE_ID,
TRT,
SEX,
bidose,
NO_TOT_DIST,
OF_TOT_DIST,
NS_TOT_DIST,
SP_TOT_DIST,
PP_TOT_DIST,
NS_DUR_STRANGER,
SP_DUR_STRANGER,
PP_STRANGER_TOT_DUR
) %>% na.omit()
BehDat <- BehDat.cntrl
N <- ncol(BehDat.cntrl)
Summary Visualizations
#Reorder Treatment Levels
BehDat.cntrl$TRT <- factor(BehDat.cntrl$TRT, levels = c("CTRL", "LOW", "MID", "HIGH"))
#Create Boxplots
for(i in 5:N) {
a <- ggplot(BehDat.cntrl,aes(x=bidose,y=BehDat.cntrl[,i])) +
geom_boxplot() +
geom_jitter(aes(color=SEX)) +
ggtitle(colnames(BehDat.cntrl)[i]) +
facet_wrap(~SEX)
print(a)
}
#Create Boxplots
for(i in 5:N) {
a <- ggplot(BehDat.cntrl,aes(x=TRT,y=BehDat.cntrl[,i])) +
geom_boxplot() +
geom_jitter(aes(color=SEX)) +
ggtitle(colnames(BehDat.cntrl)[i]) +
facet_wrap(~SEX)
print(a)
}
#Create Histograms
for(i in 5:N) {
a <- ggplot(BehDat.cntrl,aes(scale(BehDat.cntrl[,i]))) +
geom_histogram(binwidth = .5) +
ggtitle(colnames(BehDat.cntrl)[i])
print(a)
}
#Create Histograms
for(i in 5:N) {
a <- ggplot(BehDat.cntrl,aes(BehDat.cntrl[,i])) +
geom_histogram() +
ggtitle(colnames(BehDat.cntrl)[i])
print(a)
}
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#QQ Plots of Univariate Normality
for(i in 5:N) {
qqnorm(scale(BehDat[,i]), main = colnames(BehDat)[i])
abline(0,1)
}
Normality Assumptions
#Univariate Normality Tests
mvn(BehDat.cntrl[,5:N])$univariateNormality
## Test Variable Statistic p value Normality
## 1 Shapiro-Wilk NO_TOT_DIST 0.9343 7e-04 NO
## 2 Shapiro-Wilk OF_TOT_DIST 0.9867 0.6126 YES
## 3 Shapiro-Wilk NS_TOT_DIST 0.9813 0.3251 YES
## 4 Shapiro-Wilk SP_TOT_DIST 0.9598 0.0167 NO
## 5 Shapiro-Wilk PP_TOT_DIST 0.9500 0.0047 NO
## 6 Shapiro-Wilk NS_DUR_STRANGER 0.8885 <0.001 NO
## 7 Shapiro-Wilk SP_DUR_STRANGER 0.8806 <0.001 NO
## 8 Shapiro-Wilk PP_STRANGER_TOT_DUR 0.9493 0.0043 NO
#Multivariate Normality Tests Full Data
mvn(BehDat.cntrl[,5:N])$multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 224.591409017898 2.37327035079921e-08 NO
## 2 Mardia Kurtosis 2.89212860767036 0.00382641279175866 NO
## 3 MVN <NA> <NA> NO
#MultiVariate Outliers
psych::outlier(BehDat.cntrl[,5:N])
## 1 2 5 6 7 8 9 10
## 4.924770 7.206021 11.313192 8.810965 11.349433 12.174015 16.084224 3.798471
## 11 12 13 15 16 17 18 19
## 24.084510 20.269713 8.539640 5.757653 4.122267 2.460637 8.218868 6.337365
## 20 21 22 23 24 25 26 28
## 5.052933 3.033441 6.556188 4.901602 5.091132 9.239003 7.122125 3.792804
## 29 31 32 33 34 35 36 37
## 9.058182 5.000235 3.321862 5.982680 13.625658 5.945709 3.530886 5.109578
## 38 39 40 42 43 44 45 46
## 6.317432 3.233111 7.777143 2.803873 6.738061 7.739162 4.144873 1.331467
## 47 48 49 50 51 53 54 55
## 6.301996 3.140668 1.594883 3.635818 5.318109 3.289004 3.190527 7.757013
## 56 57 58 60 61 62 63 65
## 9.350242 2.967024 14.972386 8.426469 5.021785 9.743765 19.651095 10.483957
## 66 67 68 69 71 72 73 74
## 5.898630 4.917003 5.256222 6.629079 6.331870 7.749944 18.707813 3.874637
## 75 76 77 78 79 80 81 82
## 10.503296 19.055410 10.194123 5.100855 18.982783 5.728210 9.269426 13.703162
## 83 84 86 87
## 16.471481 4.476908 10.423443 9.980082
md <- mahalanobis(BehDat.cntrl[,5:N], center = colMeans(BehDat.cntrl[,5:N]), cov = cov(BehDat.cntrl[,5:N]))
alpha <- .05
cutoff <- (qchisq(p = 1 - alpha, df = ncol(BehDat.cntrl[,5:N])))
names_outliers_MH <- which(md > cutoff)
excluded_mh <- names_outliers_MH
BehDat.clean <- BehDat.cntrl[-excluded_mh, ]
BehDat.cntrl[excluded_mh, ]
## VOLE_ID TRT SEX bidose NO_TOT_DIST OF_TOT_DIST NS_TOT_DIST SP_TOT_DIST
## 9 100 LOW F TRT 14586.32 89651.41 11403.55 6357.65
## 11 240 LOW M TRT 13247.24 50260.21 7017.38 7776.55
## 12 115 MID M TRT 9206.03 66547.37 12071.38 5397.70
## 63 118 LOW M TRT 5079.34 84551.70 970.11 6150.77
## 73 25 CTRL F CTRL 34146.01 99487.90 7711.84 13615.56
## 76 239 CTRL F CTRL 37148.55 102330.20 12145.68 17682.69
## 79 140 HIGH F TRT 5318.38 108014.24 6476.99 6994.57
## 83 197 HIGH M TRT 6493.06 42034.91 3527.24 2764.85
## PP_TOT_DIST NS_DUR_STRANGER SP_DUR_STRANGER PP_STRANGER_TOT_DUR
## 9 24046.73 133.53 171.14 186.72
## 11 12960.20 204.64 436.27 202.54
## 12 10750.26 29.53 122.86 416.58
## 63 19812.46 190.52 102.34 204.50
## 73 10489.75 33.90 65.53 212.21
## 76 16740.04 1.77 82.75 120.25
## 79 7938.34 35.90 94.63 46.25
## 83 7867.31 414.78 166.20 52.42
psych::outlier(BehDat.clean[,5:N])
## 1 2 5 6 7 8 10 13
## 6.634019 9.483911 13.785646 10.831975 14.614310 14.416972 3.882590 9.415907
## 15 16 17 18 19 20 21 22
## 7.842837 5.080425 2.436416 10.521681 6.296086 5.163644 3.021169 6.869630
## 23 24 25 26 28 29 31 32
## 5.495161 5.865175 11.884975 6.965591 3.813394 12.696770 5.120620 4.169479
## 33 34 35 36 37 38 39 40
## 8.916738 18.059859 6.641050 6.023411 5.836637 6.979330 3.510208 10.171394
## 42 43 44 45 46 47 48 49
## 2.922914 6.632085 8.470148 5.505343 1.613534 6.712723 3.196198 2.358368
## 50 51 53 54 55 56 57 58
## 5.089751 5.093937 4.811738 3.742961 7.644894 13.102781 4.649964 15.253012
## 60 61 62 65 66 67 68 69
## 10.883619 5.540172 9.678966 11.477874 7.222190 5.306813 8.468590 11.628350
## 71 72 74 75 77 78 80 81
## 7.956107 8.797285 5.030881 11.921585 14.355710 6.326730 6.805216 12.072400
## 82 84 86 87
## 14.977312 6.619728 11.113004 10.570112
#Univariate Normality Tests w/o Outliers
mvn(BehDat.clean[,5:N])$univariateNormality
## Test Variable Statistic p value Normality
## 1 Shapiro-Wilk NO_TOT_DIST 0.9720 0.1296 YES
## 2 Shapiro-Wilk OF_TOT_DIST 0.9864 0.6684 YES
## 3 Shapiro-Wilk NS_TOT_DIST 0.9701 0.1013 YES
## 4 Shapiro-Wilk SP_TOT_DIST 0.9666 0.0651 YES
## 5 Shapiro-Wilk PP_TOT_DIST 0.9519 0.0106 NO
## 6 Shapiro-Wilk NS_DUR_STRANGER 0.8881 <0.001 NO
## 7 Shapiro-Wilk SP_DUR_STRANGER 0.8834 <0.001 NO
## 8 Shapiro-Wilk PP_STRANGER_TOT_DUR 0.9522 0.0109 NO
#Multivariate Tests W/o Outliers
mvn(BehDat.clean[,5:N])$multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 125.236735694877 0.353410002939259 YES
## 2 Mardia Kurtosis -0.626276555952885 0.531133560791075 YES
## 3 MVN <NA> <NA> YES
Correlations
#Full Data
res2.A <-rcorr(as.matrix(BehDat.cntrl[,5:N]),type = "pearson")
corrplot(res2.A$r, type="upper", order="hclust",tl.cex=.8)
#Clean Data
res2.A <-rcorr(as.matrix(BehDat.clean[,5:N]),type = "pearson")
corrplot(res2.A$r, type="upper", order="hclust")
#Full Data Male
res2.A <-rcorr(as.matrix(BehDat.cntrl[BehDat.cntrl$SEX == 'M',5:N]),type = "pearson")
corrplot(res2.A$r, type="upper", order="hclust", title = 'Full Data Male')
#Full Data Female
res2.A <-rcorr(as.matrix(BehDat.cntrl[BehDat.cntrl$SEX == 'F',5:N]),type = "pearson")
corrplot(res2.A$r, type="upper", order="hclust", title = 'Full Data Female')
#Table Count
table(BehDat.cntrl$bidose,BehDat.cntrl$SEX)
##
## F M
## CTRL 10 11
## TRT 28 27
#Correlations of Sex and Exposure
res.1 <-rcorr(as.matrix(filter(BehDat.cntrl, SEX == 'M' & bidose == 'CTRL')[,5:N]),type = "pearson")
res.2 <-rcorr(as.matrix(filter(BehDat.cntrl, SEX == 'M' & bidose == 'TRT')[,5:N]),type = "pearson")
res.3 <-rcorr(as.matrix(filter(BehDat.cntrl, SEX == 'F' & bidose == 'CTRL')[,5:N]),type = "pearson")
res.4 <-rcorr(as.matrix(filter(BehDat.cntrl, SEX == 'F' & bidose == 'TRT')[,5:N]),type = "pearson")
#Corrplots of Sex and Exposure
corrplot(res.1$r, type="upper" , title = 'MALE CTRL',tl.cex=.8,mar=c(0,0,2,0))
corrplot(res.2$r, type="upper" , title = 'MALE TRT',tl.cex=.8,mar=c(0,0,2,0))
corrplot(res.3$r, type="upper" , title = 'FEMALE CTRL',tl.cex=.8,mar=c(0,0,2,0))
corrplot(res.4$r, type="upper" , title = 'FEMALE TRT',tl.cex=.8,mar=c(0,0,2,0))
Exploratory Factor Analysis
#Inputs
#Preprocessing
X <- BehDat.cntrl[,5:N]
X.c <- BehDat.clean[,5:N]
for (i in 1:5) {
X[,i] <- max(X[,i])-X[,i]
X.c[,i] <- max(X.c[,i])-X.c[,i]
}
X <- scale(X)
X.c <- scale(X.c)
#Parallel Analysis
fa.parallel(X, n.iter =100, fm = "ml", fa="fa")
## Parallel analysis suggests that the number of factors = 2 and the number of components = NA
fa.parallel(X.c, n.iter =100, fm = "ml", fa="fa")
## Parallel analysis suggests that the number of factors = 2 and the number of components = NA
f <- 2
#EFA Full Data
fa.out.none <- fa(X, fm="ml",nfactors = f, rotate="none")
fa.out.varimax <- fa(X, fm="ml",nfactors = f, rotate="varimax")
fa.out.quartimax <- fa(X, fm="ml",nfactors = f, rotate="quartimax")
Results <- rbind(fa.out.none$TLI,fa.out.none$PVAL,fa.out.none$RMSEA[1])
rownames(Results) <- c("Tucker Lewis Index", "ChiSq Pval","RMSEA")
colnames(Results) <- "Fit"
round(Results,2)
## Fit
## Tucker Lewis Index 1.06
## ChiSq Pval 0.87
## RMSEA 0.00
par(mfrow= c(1,3))
fa.diagram(fa.out.none, cut = 0.2, simple = F, main = "No Rotation")
fa.diagram(fa.out.varimax, cut =0.2, simple = F, main = "Varimax rotation")
fa.diagram(fa.out.quartimax,cut =0.16, simple = F, main = "Qaurtimax roation")
#EFA Fo Outliers
fa.out.none <- fa(X.c, fm="ml",nfactors = f, rotate="none")
fa.out.varimax <- fa(X.c, fm="ml",nfactors = f, rotate="varimax")
fa.out.quartimax <- fa(X.c, fm="ml",nfactors = f, rotate="quartimax")
Results <- rbind(fa.out.none$TLI,fa.out.none$PVAL,fa.out.none$RMSEA[1])
rownames(Results) <- c("Tucker Lewis Index", "ChiSq Pval","RMSEA")
round(Results,2)
## RMSEA
## Tucker Lewis Index 0.92
## ChiSq Pval 0.08
## RMSEA 0.09
par(mfrow= c(1,3))
fa.diagram(fa.out.none, cut = 0.2, simple = F, main = "No Rotation")
fa.diagram(fa.out.varimax, cut =0.2, simple = F, main = "Varimax rotation")
fa.diagram(fa.out.quartimax, cut=0.2, simple = F, main = "Qaurtimax roation")
Confirmatory Factor Analysis (Total Stranger Interactions)
#Preprocessing
CFADat <- BehDat.cntrl[,5:N]
CFADat.c <- BehDat.clean[,5:N]
for (i in 1:5) {
CFADat[,i] <- max(CFADat[,i])-CFADat[,i]
CFADat.c[,i] <- max(CFADat.c[,i])-CFADat.c[,i]
}
CFADat <- data.frame(scale(CFADat))
CFADat.c <- data.frame(scale(CFADat.c))
#Model Specification
textModel <- sem::specifyModel(text="
##Specification of Anxiety
ANX -> NO_TOT_DIST, lambda11, NA
ANX -> OF_TOT_DIST, lambda21, NA
ANX -> NS_TOT_DIST, lambda31, NA
ANX -> PP_TOT_DIST, lambda41, NA
ANX -> SP_TOT_DIST, lambda51, NA
##Specification of Social
SOC -> NS_DUR_STRANGER, lambda12, NA
SOC -> PP_STRANGER_TOT_DUR, lambda22, NA
SOC -> SP_DUR_STRANGER, lambda32, NA
##Uniqueness of Variables
NO_TOT_DIST <-> NO_TOT_DIST, psi1, NA
OF_TOT_DIST <-> OF_TOT_DIST, psi2, NA
NS_TOT_DIST <-> NS_TOT_DIST, psi3, NA
PP_TOT_DIST <-> PP_TOT_DIST, psi4, NA
SP_TOT_DIST <-> SP_TOT_DIST, psi5, NA
NS_DUR_STRANGER <-> NS_DUR_STRANGER, psi6, NA
PP_STRANGER_TOT_DUR <-> PP_STRANGER_TOT_DUR, psi7, NA
SP_DUR_STRANGER <-> SP_DUR_STRANGER, psi8, NA
##Assumed Fixed Variances of Factors (since data is scaled?)
ANX <-> ANX, NA, 1
SOC <-> SOC, NA, 1
##Correlation Between Latent Variables
ANX <-> SOC, rho1, NA
")
## NOTE: it is generally simpler to use specifyEquations() or cfa()
## see ?specifyEquations
#CFA for All Data
Risk_Model1 <- sem::sem(model = textModel,
data = CFADat,
standardized = T,
objective=objectiveML
)
# Model Chisquare = 22.09628 Df = 19 Pr(>Chisq) = 0.2794909
# Goodness-of-fit index = 0.9334547
# Adjusted goodness-of-fit index = 0.8739142
# RMSEA index = 0.04661358 90% CI: (NA, 0.1155219)
# Bentler-Bonett NFI = 0.9031891
# Tucker-Lewis NNFI = 0.9772129
#1 Estimate insig (8.763709e-02 SP_DUR_STRANGER <--> SP_DUR_STRANGER)
#CFA for All Data
Risk_Model2 <- sem::sem(model = textModel,
data = CFADat,
standardized = T,
objective=objectiveGLS
)
# Model Chisquare = 18.69105 Df = 19 Pr(>Chisq) = 0.4768105
# Goodness-of-fit index = 0.9376965
# Adjusted goodness-of-fit index = 0.8819513
# RMSEA index = 0 90% CI: (NA, 0.09925957)
# Bentler-Bonett NFI = 0.9697971
# Tucker-Lewis NNFI = 1.000771
#No insig Estimates
#CFA w/o outliers
Risk_Model3 <- sem::sem(model = textModel,
data = CFADat.c,
standardized = T,
objective=objectiveML
)
# Model Chisquare = 41.12428 Df = 19 Pr(>Chisq) = 0.002323121
# Goodness-of-fit index = 0.8674984
# Adjusted goodness-of-fit index = 0.7489442
# RMSEA index = 0.1318319 90% CI: (0.07607498, 0.187205)
# Bentler-Bonett NFI = 0.8451127
# Tucker-Lewis NNFI = 0.8627255
#CFA w/o outliers
Risk_Model4 <- sem::sem(model = textModel,
data = CFADat.c,
standardized = T,
objective=objectiveGLS
)
# Model Chisquare = 28.62445 Df = 19 Pr(>Chisq) = 0.07212646
# Goodness-of-fit index = 0.8931924
# Adjusted goodness-of-fit index = 0.7976276
# RMSEA index = 0.08695082 90% CI: (NA, 0.1486247)
# Bentler-Bonett NFI = 0.9793594
# Tucker-Lewis NNFI = 0.9895619
#Fit and Summary
summary(Risk_Model2, fit.indices=c("GFI", "AGFI", "RMSEA", "NFI", "NNFI", "CFI", "RNI",
"IFI", "SRMR", "AIC", "AICc", "BIC", "CAIC"))
##
## Model Chisquare = 18.69105 Df = 19 Pr(>Chisq) = 0.4768105
## Goodness-of-fit index = 0.9376965
## Adjusted goodness-of-fit index = 0.8819513
## RMSEA index = 0 90% CI: (NA, 0.09925957)
## Bentler-Bonett NFI = 0.9697971
## Tucker-Lewis NNFI = 1.000771
## Bentler CFI = 1
## Bentler RNI = 1.000523
## Bollen IFI = 1.000515
## SRMR = 0.08891685
##
## Normalized Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.99434 -0.12222 -0.02868 0.18420 0.35386 2.11986
##
## R-square for Endogenous Variables
## NO_TOT_DIST OF_TOT_DIST NS_TOT_DIST PP_TOT_DIST
## 0.4990 0.4069 0.5565 0.5701
## SP_TOT_DIST NS_DUR_STRANGER PP_STRANGER_TOT_DUR SP_DUR_STRANGER
## 0.7795 0.5738 0.0848 0.6691
##
## Parameter Estimates
## Estimate Std Error z value Pr(>|z|)
## lambda11 0.6524931 0.12765573 5.111350 3.198647e-07
## lambda21 0.6197318 0.11608358 5.338669 9.363161e-08
## lambda31 0.7329260 0.10913082 6.716031 1.867412e-11
## lambda41 0.7397119 0.10730244 6.893709 5.435599e-12
## lambda51 0.8730125 0.09965670 8.760199 1.949066e-18
## lambda12 0.6686076 0.13826382 4.835738 1.326524e-06
## lambda22 0.2849193 0.13552158 2.102391 3.551905e-02
## lambda32 0.7446642 0.13350053 5.577987 2.433178e-08
## psi1 0.4274297 0.09176158 4.658047 3.192240e-06
## psi2 0.5598372 0.10412107 5.376791 7.582506e-08
## psi3 0.4280626 0.08882426 4.819208 1.441294e-06
## psi4 0.4126696 0.08365710 4.932870 8.103014e-07
## psi5 0.2156090 0.06779553 3.180283 1.471311e-03
## psi6 0.3320531 0.11658089 2.848263 4.395855e-03
## psi7 0.8758049 0.15076062 5.809242 6.275642e-09
## psi8 0.2742979 0.12886074 2.128638 3.328419e-02
## rho1 -0.6182661 0.13386833 -4.618464 3.865909e-06
##
## lambda11 NO_TOT_DIST <--- ANX
## lambda21 OF_TOT_DIST <--- ANX
## lambda31 NS_TOT_DIST <--- ANX
## lambda41 PP_TOT_DIST <--- ANX
## lambda51 SP_TOT_DIST <--- ANX
## lambda12 NS_DUR_STRANGER <--- SOC
## lambda22 PP_STRANGER_TOT_DUR <--- SOC
## lambda32 SP_DUR_STRANGER <--- SOC
## psi1 NO_TOT_DIST <--> NO_TOT_DIST
## psi2 OF_TOT_DIST <--> OF_TOT_DIST
## psi3 NS_TOT_DIST <--> NS_TOT_DIST
## psi4 PP_TOT_DIST <--> PP_TOT_DIST
## psi5 SP_TOT_DIST <--> SP_TOT_DIST
## psi6 NS_DUR_STRANGER <--> NS_DUR_STRANGER
## psi7 PP_STRANGER_TOT_DUR <--> PP_STRANGER_TOT_DUR
## psi8 SP_DUR_STRANGER <--> SP_DUR_STRANGER
## rho1 SOC <--> ANX
##
## Iterations = 20
round(Risk_Model2$S - Risk_Model2$C,2)
## NO_TOT_DIST OF_TOT_DIST NS_TOT_DIST SP_TOT_DIST PP_TOT_DIST
## NO_TOT_DIST 0.15 0.04 -0.02 0.02 0.03
## OF_TOT_DIST 0.04 0.06 -0.01 -0.01 0.05
## NS_TOT_DIST -0.02 -0.01 0.03 -0.02 -0.01
## SP_TOT_DIST 0.02 -0.01 -0.02 0.02 -0.01
## PP_TOT_DIST 0.03 0.05 -0.01 -0.01 0.04
## NS_DUR_STRANGER 0.24 0.05 -0.10 0.04 -0.01
## SP_DUR_STRANGER 0.19 0.13 0.00 -0.04 -0.01
## PP_STRANGER_TOT_DUR -0.03 -0.05 0.00 -0.05 -0.12
## NS_DUR_STRANGER SP_DUR_STRANGER PP_STRANGER_TOT_DUR
## NO_TOT_DIST 0.24 0.19 -0.03
## OF_TOT_DIST 0.05 0.13 -0.05
## NS_TOT_DIST -0.10 0.00 0.00
## SP_TOT_DIST 0.04 -0.04 -0.05
## PP_TOT_DIST -0.01 -0.01 -0.12
## NS_DUR_STRANGER 0.22 0.14 -0.01
## SP_DUR_STRANGER 0.14 0.17 -0.01
## PP_STRANGER_TOT_DUR -0.01 -0.01 0.04
Vhat <- Risk_Model2$vcov
corrplot(cov2cor(Vhat))
pathDiagram(Risk_Model2,
ignore.double = F,
edge.labels = "both",
file = "Risk_Model_plot",
output.type = "dot",
node.colors = c("steelblue","transparent"))
grViz("Risk_Model_plot.dot")
Factor Score Analysis
FS <- fscores(Risk_Model2)
newdat <- cbind(BehDat.cntrl[,1:4],FS)
N <- ncol(newdat)
#Reorder Treatment Levels
newdat$TRT <- factor(newdat$TRT, levels = c("CTRL", "LOW", "MID", "HIGH"))
#Create Boxplots
for(i in 5:N) {
a <- ggplot(newdat,aes(x=bidose,y=newdat[,i])) +
geom_boxplot() +
geom_jitter(aes(color=SEX)) +
ggtitle(colnames(newdat)[i]) +
facet_wrap(~SEX)
print(a)
}
#Create Boxplots
for(i in 5:N) {
a <- ggplot(newdat,aes(x=TRT,y=newdat[,i])) +
geom_boxplot() +
geom_jitter(aes(color=SEX)) +
ggtitle(colnames(newdat)[i]) +
facet_wrap(~SEX)
print(a)
}
fit1.A <- lm(ANX~TRT*SEX, data=newdat)
summary(fit1.A)
##
## Call:
## lm(formula = ANX ~ TRT * SEX, data = newdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.66054 -0.54952 0.09291 0.67205 1.72538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8586 0.2842 -3.021 0.00355 **
## TRTLOW 0.9920 0.3848 2.578 0.01211 *
## TRTMID 1.2014 0.4263 2.818 0.00632 **
## TRTHIGH 1.1609 0.4263 2.723 0.00821 **
## SEXM 0.4653 0.3927 1.185 0.24020
## TRTLOW:SEXM -0.3018 0.5679 -0.531 0.59681
## TRTMID:SEXM -0.7389 0.6243 -1.183 0.24074
## TRTHIGH:SEXM -0.5092 0.5633 -0.904 0.36923
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8987 on 68 degrees of freedom
## Multiple R-squared: 0.1861, Adjusted R-squared: 0.1023
## F-statistic: 2.222 on 7 and 68 DF, p-value: 0.04289
summary.aov(fit1.A)
## Df Sum Sq Mean Sq F value Pr(>F)
## TRT 3 11.03 3.677 4.553 0.00575 **
## SEX 1 0.23 0.232 0.287 0.59391
## TRT:SEX 3 1.30 0.432 0.535 0.65981
## Residuals 68 54.92 0.808
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot(fit1.A)
fit1.S <- lm(SOC~TRT*SEX, data=newdat)
summary(fit1.S)
##
## Call:
## lm(formula = SOC ~ TRT * SEX, data = newdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5410 -0.6372 -0.2197 0.6866 2.1874
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3891 0.2995 1.299 0.1983
## TRTLOW -0.3475 0.4055 -0.857 0.3945
## TRTMID -0.5964 0.4493 -1.328 0.1888
## TRTHIGH -1.1295 0.4493 -2.514 0.0143 *
## SEXM -0.1053 0.4138 -0.255 0.7998
## TRTLOW:SEXM 0.5062 0.5985 0.846 0.4007
## TRTMID:SEXM 0.1594 0.6580 0.242 0.8093
## TRTHIGH:SEXM 0.6499 0.5936 1.095 0.2775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9471 on 68 degrees of freedom
## Multiple R-squared: 0.133, Adjusted R-squared: 0.04371
## F-statistic: 1.49 on 7 and 68 DF, p-value: 0.1857
summary.aov(fit1.S)
## Df Sum Sq Mean Sq F value Pr(>F)
## TRT 3 7.05 2.3490 2.619 0.0579 .
## SEX 1 0.98 0.9825 1.095 0.2990
## TRT:SEX 3 1.33 0.4417 0.492 0.6887
## Residuals 68 61.00 0.8971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot(fit1.S)
fit2.A <- lm(ANX~bidose*SEX, data=newdat)
summary(fit2.A)
##
## Call:
## lm(formula = ANX ~ bidose * SEX, data = newdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.62990 -0.53938 0.09624 0.65204 1.72538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8586 0.2773 -3.096 0.00280 **
## bidoseTRT 1.1001 0.3231 3.405 0.00108 **
## SEXM 0.4653 0.3832 1.214 0.22864
## bidoseTRT:SEXM -0.4790 0.4503 -1.064 0.29102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.877 on 72 degrees of freedom
## Multiple R-squared: 0.1794, Adjusted R-squared: 0.1452
## F-statistic: 5.247 on 3 and 72 DF, p-value: 0.002484
summary.aov(fit2.A)
## Df Sum Sq Mean Sq F value Pr(>F)
## bidose 1 10.97 10.971 14.264 0.000324 ***
## SEX 1 0.27 0.266 0.346 0.558110
## bidose:SEX 1 0.87 0.870 1.131 0.291017
## Residuals 72 55.38 0.769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot(fit2.A)
fit2.S <- lm(SOC~bidose*SEX, data=newdat)
summary(fit2.S)
##
## Call:
## lm(formula = SOC ~ bidose * SEX, data = newdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5410 -0.6832 -0.2149 0.6529 2.5099
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3891 0.3031 1.283 0.2035
## bidoseTRT -0.6421 0.3531 -1.818 0.0732 .
## SEXM -0.1053 0.4188 -0.252 0.8021
## bidoseTRT:SEXM 0.3610 0.4922 0.733 0.4657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9586 on 72 degrees of freedom
## Multiple R-squared: 0.05958, Adjusted R-squared: 0.02039
## F-statistic: 1.52 on 3 and 72 DF, p-value: 0.2165
summary.aov(fit2.S)
## Df Sum Sq Mean Sq F value Pr(>F)
## bidose 1 3.23 3.235 3.520 0.0647 .
## SEX 1 0.46 0.462 0.503 0.4804
## bidose:SEX 1 0.49 0.494 0.538 0.4657
## Residuals 72 66.16 0.919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stats::anova(fit2.A,fit1.A)
## Analysis of Variance Table
##
## Model 1: ANX ~ bidose * SEX
## Model 2: ANX ~ TRT * SEX
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 72 55.376
## 2 68 54.923 4 0.4531 0.1402 0.9667
stats::anova(fit2.S,fit1.S)
## Analysis of Variance Table
##
## Model 1: SOC ~ bidose * SEX
## Model 2: SOC ~ TRT * SEX
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 72 66.164
## 2 68 61.001 4 5.1631 1.4389 0.2307
fit <- manova(cbind(ANX,SOC)~TRT*SEX, data=newdat)
summary(fit)
## Df Pillai approx F num Df den Df Pr(>F)
## TRT 3 0.254038 3.2980 6 136 0.004623 **
## SEX 1 0.050393 1.7777 2 67 0.176899
## TRT:SEX 3 0.053596 0.6242 6 136 0.710706
## Residuals 68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit <- manova(cbind(ANX,SOC)~bidose*SEX, data=newdat)
summary(fit)
## Df Pillai approx F num Df den Df Pr(>F)
## bidose 1 0.168843 7.2115 2 71 0.001408 **
## SEX 1 0.030033 1.0992 2 71 0.338746
## bidose:SEX 1 0.015586 0.5621 2 71 0.572539
## Residuals 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(fit)
## Response ANX :
## Df Sum Sq Mean Sq F value Pr(>F)
## bidose 1 10.971 10.9710 14.2644 0.0003242 ***
## SEX 1 0.266 0.2663 0.3462 0.5581100
## bidose:SEX 1 0.870 0.8702 1.1315 0.2910166
## Residuals 72 55.376 0.7691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response SOC :
## Df Sum Sq Mean Sq F value Pr(>F)
## bidose 1 3.235 3.2349 3.5202 0.06468 .
## SEX 1 0.462 0.4623 0.5031 0.48043
## bidose:SEX 1 0.494 0.4943 0.5379 0.46567
## Residuals 72 66.164 0.9189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(newdat$ANX)
hist(newdat$SOC)
